R包circlize绘制弦状图示例
鲤小白：当初学习circlize包，最大的原因是白鱼同学一直以来都很讨厌perl……
路人：等等，这俩有什么关系呢？
鲤小白：因为要画圈图啊，但是circos是perl程序，不想用，所以得找个替代品，正好看到了circlize。
路人：好吧你赢了……
白鱼同学和perl之间的梗，哪天有兴趣了可以当笑话给大家分享，足以写满一页纸的笑话了。不过这里得回过话题来，circlize包绘制弦状图。关于circlize的基础本篇不涉及，还请自行了解，本篇以一个也不算很复杂的弦状图为例，展示一下常规做法，circlize绘制的弦状图还是很好看滴。
circlize文档：https://jokergoo.github.io/circlize_book/book/
本篇所用的作图示例文件、R脚本等的百度盘链接：
https://pan.baidu.com/s/10BCdWESl6Xt7DliUyT6drw
网盘附件中的数据为某16S高通量测序所得的细菌OTU丰度、物种分类等数据。本篇将使用这些数据绘制弦状图，用以展示20种重要的OTU在6个样本中的丰度关系。20种OTU可划分为5个门（物种分类单元，界门纲目科属种的“门”），6个样本可划分为2个分组。
左侧为circlize做图结果，分为5圈：
第一圈，各OTU的门水平分类以及各样本的分组信息；
第二圈，OTU相对丰度的百分比信息；
第三圈，OTU及样本主区块，以不同颜色和标签区分，区块外周的刻度为OTU的绝对丰度信息；
第四圈，OTU及样本副区块，与主区块（第三圈）对应，展示了各OTU在各样本中的丰度，以及各样本所含每种OTU的丰度信息；
第五圈，与OTU及样本副区块（第四圈）相对应，连线展示OTU、样本关联信息。
右侧为图例，展示了各OTU（本示例中共计20种OTU）的详细物种分类详情。
准备文件
网盘附件3个文件。
（1）物种丰度表格文件（otu_table.txt）。每一列为样本，每一行为物种（此处为OTU），交叉区为个OTU在各样本中的丰度。
（2）样本分组信息文件（group.txt）。第一列为样本ID，第二列为样本分组信息。
（3）物种分类信息文件（taxonomy.txt）。第一列为OTU ID，第二列为OTU的“门水平”分类，第三列为OTU的详细物种分类。
circlize作图
预处理数据
首先根据OTU分类文件（taxonomy.txt）以及样本分组文件（group.txt）中的OTU、样本顺序，对OTU丰度表（otu_table.txt）中的OTU、样本进行排序，以及获取OTU和样本的名称。
#读取 taxonomy.txt 的内容，获取“OTU/分类”排序，OTU 名称
taxonomy <- read.delim('taxonomy.txt', sep = '\t', stringsAsFactors = FALSE)
tax_phylum <- unique(taxonomy$phylum)
taxonomy$phylum <- factor(taxonomy$phylum, levels = tax_phylum)
all_otu <- taxonomy$OTU_ID
taxonomy$OTU_ID <- factor(taxonomy$OTU_ID, levels = all_otu)
#读取 group.txt 的内容，获取“样本/分组”排序，样本名称
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)
all_group <- unique(group$group_ID)
group$group_ID <- factor(group$group_ID, levels = all_group)
all_sample <- group$sample_ID
#读取 otu_table.txt，排序 OTU 和样本
otu_table <- read.delim('otu_table.txt', sep = '\t')
otu_table <- merge(taxonomy, otu_table, by = 'OTU_ID')
otu_table <- otu_table[order(otu_table$phylum, otu_table$OTU_ID), ]
rownames(otu_table) <- otu_table$OTU_ID
otu_table <- otu_table[all_sample]
然后作些小统计，各OTU的总丰度，各样本中所有OTU的总丰度，各OTU在各样本中的丰度，以及各样本所含各OTU的丰度。顺便给各OTU及样本定义颜色。这些统计好的结果将直接用于circlize作图。
library(reshape2)##生成作图数据
#circlize 外圈属性数据
all_ID <- c(all_otu, all_sample)
accum_otu <- rowSums(otu_table)
accum_sample <- colSums(otu_table)
all_ID_xlim <- cbind(rep(0, length(all_ID)),data.frame(c(accum_otu, accum_sample)))
#circlize 内圈连线数据
otu_table$otu_ID <- all_otu
plot_data <- melt(otu_table, id = 'otu_ID') #此处使用了reshape2包中的melt()命令
colnames(plot_data)[2] <- 'sample_ID'
plot_data$otu_ID <- factor(plot_data$otu_ID, levels = all_otu)
plot_data$sample_ID <- factor(plot_data$sample_ID, levels = all_sample)
plot_data <- plot_data[order(plot_data$otu_ID, plot_data$sample_ID), ]
plot_data <- plot_data[c(2, 1, 3, 3)]
#颜色设置
color_otu <- c('#8DD3C7', '#FFFFB3', '#BEBADA', '#FB8072', '#80B1D3', '#FDB462', '#B3DE69', '#FCCDE5', '#BC80BD', '#CCEBC5', '#FFED6F', '#E41A1C', '#377EB8', '#4DAF4A', '#984EA3', '#FF7F00', '#FFFF33', '#A65628', '#F781BF', '#66C2A5')
color_sample <- c('#6181BD', '#F34800', '#64A10E', '#FF00FF', '#c7475b', '#049a0b')
color_phylum <- c('#BEAED4', '#FDC086', '#FFFF99', '#386CB0', '#F0027F')
color_group <- c('#4253ff', '#ff4308')
names(color_otu) <- all_otu
names(color_sample) <- all_sample
定义整体布局
然后接下来就是使用circlize包绘制弦状图了。
注：绘图细节部分的调整很多也很繁琐，此处默认使用的各细节参数设置（如字体大小等）均根据示例文件而来。若是大家后续更换为自己的数据时，还需多加调试参数了。
首先定义整体的布局，这一步生成空白画板。
library(circlize)#整体画板布局
gap_size <- c(rep(3, length(all_otu) - 1), 6, rep(3, length(all_sample) - 1), 6)
circos.par(cell.padding = c(0, 0, 0, 0), start.degree = 270, gap.degree = gap_size)
circos.initialize(factors = factor(all_ID, levels = all_ID), xlim = all_ID_xlim)
gap_size将用于设定circlize图中，各OTU或样本区块之间的间距（如下图所示）；circos.par就是画板设置了；circos.initialize定义绘图因子，此处使用了前述得到的circlize外圈属性数据框all_ID_xlim。
第一圈，OTU分类、样本分组区块
然后开始绘图。R的circlize也是由外圈往内圈逐个添加。
#绘制 OTU 分类、样本分组区块（第一圈）
circos.trackPlotRegion(
ylim = c(0, 1), track.height = 0.03, bg.border = NA,
panel.fun = function(x, y) {
sector.index = get.cell.meta.data('sector.index')
xlim = get.cell.meta.data('xlim')
ylim = get.cell.meta.data('ylim')
} )
for (i in 1:length(tax_phylum)) {
tax_OTU <- {subset(taxonomy, phylum == tax_phylum[i])}$OTU_ID
highlight.sector(tax_OTU, track.index = 1, col = color_phylum[i], text = tax_phylum[i], cex = 0.5, text.col = 'black', niceFacing = FALSE)
}
for (i in 1:length(all_group)) {
group_sample <- {subset(group, group_ID == all_group[i])}$sample_ID
highlight.sector(group_sample, track.index = 1, col = color_group[i], text = all_group[i], cex = 0.7, text.col = 'black', niceFacing = FALSE)
}
第一条命令用于绘制第一圈。ylim控制y轴范围，track.height控制高度，bg.border控制边框颜色；定义的函数panel.fun()以循环的方式，读取all_ID_xlim中的内容，逐一绘制外圈区块。运行完上述的第一条命令之后，发现绘图区域啥也没显示。但事实上，并非没有绘制任何元素，只是将绘制元素的填充颜色及边框颜色设置为NA了，因为此处不需要任何颜色作为添加。
然后接下来使用两个循环，将预先设定的OTU门分类颜色以及样本分组颜色添加在刚才未定义颜色的第一圈。highlight.sector()实现这个功能，并可以将处于同一分组的样本合并在一起。track.index = 1指定了将颜色添加在circlize图的第一圈；col和text参数分别指定颜色和标签内容；cex指定标签字体大小，text.col指定标签字体颜色，niceFacing用于展示标签字体的方向，默认为TRUE但个人推荐使用FALSE。更多参数信息可?highlight.sector。
第二圈，OTU丰度的百分比信息
继续往内圈绘制。
#添加 OTU 百分比注释（第二圈）circos.trackPlotRegion(
ylim = c(0, 1), track.height = 0.05, bg.border = NA,
panel.fun = function(x, y) {
sector.index = get.cell.meta.data('sector.index')
xlim = get.cell.meta.data('xlim')
ylim = get.cell.meta.data('ylim')
} )
circos.track(
track.index = 2, bg.border = NA,
panel.fun = function(x, y) {
xlim = get.cell.meta.data('xlim')
ylim = get.cell.meta.data('ylim')
sector.name = get.cell.meta.data('sector.index')
xplot = get.cell.meta.data('xplot')
by = ifelse(abs(xplot[2] - xplot[1]) > 30, 0.25, 1)
for (p in c(0, seq(by, 1, by = by))) circos.text(p*(xlim[2] - xlim[1]) + xlim[1], mean(ylim) + 0.4, paste0(p*100, '%'), cex = 0.4, adj = c(0.5, 0), niceFacing = FALSE)
circos.lines(xlim, c(mean(ylim), mean(ylim)), lty = 3)
} )
首先我们使用circos.trackPlotRegion()绘制一圈无边框、无颜色的空白区域，目的为添加的刻度标签等内容留出展示的空间。
然后使用circos.track()添加百分比刻度标签。其中，track.index = 2意为将刻度标签展示在circlize图的第二圈；定义的函数panel.fun()使用循环的方式，读取all_ID_xlim中的内容（丰度信息）并在读取后进行计算得到百分比数值，同时额外使用if判断OTU丰度的绝对数值，若数值足够大则考虑将展示更多的刻度（以25%为一刻度），若数值过小则仅展示0和100%这两个刻度；circos.lines()添加了外周的点状虚线。
第三、四圈， OTU及样本名称区块
其中第三圈为主区块，第四圈为副区块。
#绘制 OTU、样本主区块（第三圈）circos.trackPlotRegion(
ylim = c(0, 1), track.height = 0.03, bg.col = c(color_otu, color_sample), bg.border = NA, track.margin = c(0, 0.01),
panel.fun = function(x, y) {
xlim = get.cell.meta.data('xlim')
sector.name = get.cell.meta.data('sector.index')
circos.axis(h = 'top', labels.cex = 0.4, major.tick.percentage = 0.4, labels.niceFacing = FALSE)
circos.text(mean(xlim), 0.2, sector.name, cex = 0.4, niceFacing = FALSE, adj = c(0.5, 0))
} )
#绘制 OTU、样本副区块（第四圈）
circos.trackPlotRegion(ylim = c(0, 1), track.height = 0.03, track.margin = c(0, 0.01))
继续使用circos.trackPlotRegion()循环读取all_ID_xlim中的内容并绘制，此时可以将预先定义的OTU颜色及样本颜色添加上。同时在内部使用circos.axis()和circos.text()函数，在循环绘图的同时，将刻度线和标签添加上。此处，sector.name为在循环过程中每次读取的OTU或样本名称；circos.text()中的mean(xlim)和0.2分别意为将给定的标签文字添加在每个区块的X轴中间位置以及Y=0.2的位置。其余参数项类似上述不再多说，更细致的参数可参阅帮助。
第三圈主圈绘制完成，第四圈副圈暂未添加颜色（将在后续操作中添加）。
最内圈，连线区
最内圈为样本-OTU丰度关联信息，以连线的方式展示。
#绘制 OTU-样本关联连线（最内圈）
for (i in seq_len(nrow(plot_data))) {
circos.link(
plot_data[i,2], c(accum_otu[plot_data[i,2]], accum_otu[plot_data[i,2]] - plot_data[i,4]),
plot_data[i,1], c(accum_sample[plot_data[i,1]], accum_sample[plot_data[i,1]] - plot_data[i,3]),
col = paste0(color_otu[plot_data[i,2]], '70'), border = NA )
circos.rect(accum_otu[plot_data[i,2]], 0, accum_otu[plot_data[i,2]] - plot_data[i,4], 1, sector.index = plot_data[i,2], col = color_sample[plot_data[i,1]], border = NA)
circos.rect(accum_sample[plot_data[i,1]], 0, accum_sample[plot_data[i,1]] - plot_data[i,3], 1, sector.index = plot_data[i,1], col = color_otu[plot_data[i,2]], border = NA)
accum_otu[plot_data[i,2]] = accum_otu[plot_data[i,2]] - plot_data[i,4]
accum_sample[plot_data[i,1]] = accum_sample[plot_data[i,1]] - plot_data[i,3]
}
此处以循环的方式添加最内圈连线，连线颜色以OTU的预设颜色为准（绘制的同时将颜色的不透明度设置为70%）。先前绘制外圈第四圈的时候未添加颜色，也在此处将颜色添加上，即命令中使用的两个circos.rect()函数分别添加。
运行完毕后，我们的circlize初图就完成了。
可选，添加图例
其实到了这一步就已经完成了，不过示例图中还包含了图例信息。circlize包中是没有绘制图例的函数的，因此还需额外调用其它的包绘制图例。
在circlize的官方说明文档中，提到可通过ComplexHeatmap包Legend()函数，为circlize图添加图例。
备注：为了更好的体验，请使用3.6版本及以上的R。经过多次踩坑，发现像3.3、3.4等版本的R中，ComplexHeatmap包的很多功能无法适用，函数残缺且不兼容等问题…….
library(ComplexHeatmap) #可用此包添加图例
library(grid) #可用此包调整画板
#添加图例
otu_legend <- Legend(
at = all_otu, labels = taxonomy$detail, labels_gp = gpar(fontsize = 8),
grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'), type = 'points', pch = NA, background = color_otu)
pushViewport(viewport(x = 0.85, y = 0.5))
grid.draw(otu_legend)
upViewport()
#推荐作图完成后，清除目前的 circlize 样式，以免影响继续作图
circos.clear()
图例添加完毕后，最终结果如下所示。
其他
如果图片比例失调，除了更改作图参数外，还可通过在本地输出为pdf/png时设置保存图形的长宽比列来解决。
pdf('circlize_plot.pdf', width = 20, height = 8)
circle_size = unit(1, 'snpc')
#然后作图
#整体布局
gap_size <- c(rep(3, length(all_otu) - 1), 6, rep(3, length(all_sample) - 1), 6)
circos.par(cell.padding = c(0, 0, 0, 0), start.degree = 270, gap.degree = gap_size)
circos.initialize(factors = factor(all_ID, levels = all_ID), xlim = all_ID_xlim)
#中间那些作图步骤（略）
……
##添加图例
otu_legend <- Legend(
at = all_otu, labels = taxonomy$detail, labels_gp = gpar(fontsize = 8),
grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'), type = 'points', pch = NA, background = color_otu)
pushViewport(viewport(x = 0.85, y = 0.5))
grid.draw(otu_legend)
upViewport()
##清除 circlize 样式并关闭画板/pdf 文件
circos.clear()
dev.off()